% this code aggregates up the frequency, size and IO data using BEA
% classification

clearvars -except wd 
cd(fullfile(wd, 'Figure4_LHS'));

%cd '/Users/raphaelschoenle/Dropbox/Networks/code/JME_post/Figure4_LHS'

% initial number of sectors
n_sectors = 341;

%% initial data
% read in aggregations 1-5
% aggregation 1
agg_1 = fopen('agg_1.csv');
fmt = '%s%s';
data_1 = textscan(agg_1,fmt,'headerlines',1,'delimiter',',');
% aggregation 2
agg_2 = fopen('agg_2.csv');
fmt = '%s%s';
data_2 = textscan(agg_2,fmt,'headerlines',1,'delimiter',',');
% aggregation 3
agg_3 = fopen('agg_3.csv');
fmt = '%s%s';
data_3 = textscan(agg_3,fmt,'headerlines',1,'delimiter',',');
% aggregation 4
agg_4 = fopen('agg_4.csv');
fmt = '%s%s';
data_4 = textscan(agg_4,fmt,'headerlines',1,'delimiter',',');
% aggregation 5
agg_5 = fopen('agg_5.csv');
fmt = '%s%s';
data_5 = textscan(agg_5,fmt,'headerlines',1,'delimiter',',');

% BEA codes: these are always 341 identical codes from the first column
bea_codes = data_1{1};
% codes agg1 etc.
agg1_codes = data_1{2};
agg2_codes = data_2{2};
agg3_codes = data_3{2};
agg4_codes = data_4{2};
agg5_codes = data_5{2};

% BEA consumption data
Ck = csvread("2002CkDetailed.csv");
Omega_C = Ck/sum(Ck); % Share of consumption; the weights omega

% Data
load FPA.mat
FPA_BEA=FPA;

% sectoral IO weights 
Znum = csvread("2002RevShareDetailed.csv");
n = length(Ck); % Number of Industries
Z = Znum(1:n,1:n);
Omega = Z./(ones(n,1)*sum(Z)); %Normalize the share of inputs for each sector, omega

% I/O matrix in terms of values
total_inputs = xlsread('intermediate_use.xlsx');

% set weights to 0 when aggregating frequencies
Omega_C_W = Omega_C;
Omega_C_W(Omega_C_W<0)=0;
Omega_C_W=Omega_C_W./sum(Omega_C_W);

% loop through all 5 levels of aggregation: from 341 (data_1) to data_k
for k = 1:5
    % kth aggregation: second column is the more aggregated calibration
    aggdata = eval(['data_' num2str(k) '{2}']);
    C = setdiff(aggdata, data_1{1});

    % positions of all aggregating codes
    loc_all = 0;
    for jj=1:length(C)
        loc_all = [loc_all; find(strcmp(aggdata,C(jj)))];
    end
    loc_all = loc_all(2:end);
    % common positions
    loc_common = setdiff(1:341,loc_all);

    % aggregate up consumption weights
    for jj=1:length(C)
        loc = find(strcmp(aggdata,C(jj)));
        if jj==1
            Omega_C_new = [ Omega_C(loc_common); sum(Omega_C(loc))];
        else
            Omega_C_new = [ Omega_C_new; sum(Omega_C(loc))];
        end
    end

    % aggregate up frequency of price changes (based on BEA)
    % make sure no negative FPAs!
    for jj=1:length(C)
        loc = find(strcmp(aggdata,C(jj)));
        if jj==1
               if  sum(Omega_C_W(loc)) == 0
                FPA_BEA_new = [ FPA_BEA(loc_common); mean(FPA_BEA(loc))   ];
               else
                FPA_BEA_new = [ FPA_BEA(loc_common); (FPA_BEA(loc)'*( Omega_C_W(loc)./sum(Omega_C_W(loc))))];
               end              
        else
            FPA_BEA_new = [ FPA_BEA_new; (FPA_BEA(loc)'*Omega_C(loc))./sum(Omega_C(loc))];
        end
    end

    % aggregate up I/O matrix
    Omega=Omega';
    % columns: simply additive
    for jj=1:length(C)
        loc = find(strcmp(aggdata,C(jj)));
        if jj==1
                Omega_new = [ Omega(:,loc_common) sum(Omega(:,loc),2) ];
        else
                Omega_new = [ Omega_new sum(Omega(:,loc),2)];
        end
    end
    loc_sum = 0;
    % rows: weighted average of rows, weights are total input use based
    for jj=1:length(C)
        loc = find(strcmp(aggdata,C(jj)));
        loc_sum = loc_sum+length(loc);
        if jj==1
            Omega_new2 = [ Omega_new(loc_common,:); (total_inputs(loc)./sum(total_inputs(loc)))'...
            *Omega_new(loc,:) ];
        else
            Omega_new2(length(loc_common)+jj,:) = (total_inputs(loc)./sum(total_inputs(loc)))'...
            *Omega_new(loc,:);
        end
    end

    Omega_new2= Omega_new2';
    
    % aggregate up total input usage
    for jj=1:length(C)
        loc = find(strcmp(aggdata,C(jj)));
        if jj==1
            total_inputs_new = [ total_inputs(loc_common); sum(total_inputs(loc))];
        else
            total_inputs_new = [ total_inputs_new; sum(total_inputs(loc))];
        end
    end
    
    % renormalize Omega matrix
    Omega_new2=Omega_new2./(ones(length(Omega_C_new),1)*sum(Omega_new2)); %Normalize the share of inputs for each sector, omega
    
    % check: is the sum of input weights = 1?
    sum(Omega_new2)
    
    % save output for code to run
	filename = [wd '/Figure4_LHS/empirical_aggregation' num2str(length(Omega_C_new)) '.mat']
    save(filename, 'Omega_new2', 'FPA_BEA_new', 'Omega_C_new', 'total_inputs_new');

end
